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Abstract 

We  discuss  several  models  for  a  flash-heat  experiment  in  homogeneous 
isotropic  media.  We  use  these  to  investigate  the  use  of  homogenization  tech¬ 
niques  in  approximating  models  for  interrogation  via  flash-heating  in  porous 
materials.  We  represent  porous  materials  as  both  randomly  perforated  do¬ 
mains  and  periodically  perforated  domains. 
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1  INTRODUCTION 


1  Introduction 

Nondestructive  evaluation  (NDE)  is  an  important  tool  in  many  fields.  Nonde¬ 
structive  evaluation  is  used  to  identify  damage  in  structures,  including  com¬ 
ponents  of  aircraft,  spacecraft,  automobiles,  trains  and  piping,  as  they  age 
beyond  their  design  life.  Proper  use  of  nondestructive  evaluation  can  increase 
the  safety  and  service  life  of  components  of  many  structures.  Ultrasound, 
magnetic  particle  imaging,  eddy  current,  acoustic  emission  and  radiology  are 
a  few  examples  of  NDE  techniques.  These  techniques  have  been  developed  in  a 
large  number  of  applications,  particularly  for  homogeneous  metallic  materials 
[1]— [5] .  Composite  materials  are  increasingly  popular  in  structures  including 
aeronautical  and  aerospace  vehicles  and  are  made  with  an  acceptable  level  of 
porosity.  This  porosity  causes  a  certain  amount  of  noise  in  the  interrogation 
signal  of  any  NDE  technique.  In  [6]  where  we  investigated  using  active  ther¬ 
mography  to  detect  damage  in  porous  domains  we  found  that  we  could  detect 
damage  of  a  certain  size  within  a  porous  medium  using  thermal  interrogation. 
However,  the  methods  developed  in  [6]  were  too  computationally  intensive  to 
use  in  the  sophisticated  parameter  estimation  routines  needed  to  characterize 
damage.  As  in  [6]  we  treat  here  the  problem  of  modeling  the  flow  of  heat  in  a 
porous  domain  but  in  the  current  work  we  focus  on  the  results  of  homogeniza¬ 
tion  theory  (which  are  less  computationally  intensive  than  the  methods  used 
in  [6])  to  capture  the  behavior. 

We  aim  to  use  homogenization  theory  to  model  a  two  dimensional  cross 
section  version  of  a  flash-heat  experiment  in  porous  domains.  In  Section  2  we 
discuss  the  flash-heat  experiment  in  the  simple  case  of  isotropic  flow  in  a  ho¬ 
mogeneous  domain  and  consider  several  models  of  the  flash-heat  experiment. 
In  particular,  we  present  here  four  different  scenarios:  (i)  heat  in  a  perfectly 
insulated  two  dimensional  sample  (system  (2)  below);  (ii)  heat  flow  in  a  cross 
section  with  part  of  the  boundary  held  at  a  fixed  temperature  (system  (4)  be¬ 
low);  (iii)  heat  flow  in  a  cross  section  with  loss  only  in  the  direction  orthogonal 
to  the  cross  section  with  boundaries  otherwise  insulated  (system  (9)  below); 
and  (iv)  heat  flow  in  a  cross  section  with  small  loss  on  either  part  or  all  of 
its  boundary  (system  (12)  and  (13),  respectively,  below).  After  choosing  two 


2  SYSTEMS 


3 


models  ((9)  and  (13))  to  investigate  in  detail,  we  discuss  the  procedures  used 
in  homogenization  to  approximate  the  flow  of  heat  over  a  domain  with  random 
pores  by  the  flow  of  heat  over  a  domain  with  many  periodically  placed  pores. 
One  passes  this  periodic  domain  to  a  limit  using  homogenization  theory  to  gain 
a  limit  partial  differential  equation  which  replaces  the  periodically  perforated 
domain  with  a  homogeneous  domain  that  has  anisotropic  flow  which  approx¬ 
imates  the  effect  of  the  porosity.  With  the  results  of  homogenization  theory 
developed  in  Section  3  we  move  on  to  summarize  the  results  of  simulations 
of  the  flash-heat  experiment  on  a  randomly  perforated  domain,  a  periodically 
perforated  domain,  and  a  homogeneous  domain  with  the  anisotropic  flow  de¬ 
rived  from  homogenization  theory  in  Section  4.  We  also  graphically  analyze 
the  two  example  simulations  of  primary  interest  in  Section  4. 


2  Systems 


We  first  recall  the  system  (1)  which  corresponds  to  the  physical  flash-heat 
experiment  described  in  [6].  This  experiment  assumes  that  the  temperature  of 
the  specimen  is  within  the  solid  state  phase  and  the  boundaries  are  perfectly 
insulated.  We  model  the  flash-heat  experiment  on  an  L-2  (length  in  the  y 
direction)  by  L\  (length  in  the  x  direction)  rectangle  during  the  time  interval 
t  e  (0,  T)  with  T  <  oo.  We  refer  to  the  L\  by  L-2  rectangle  as  and  the 
four  boundaries  u y,  for  i  G  {1,2,  3, 4}.  When  referring  to  the  entire  boundary 
of  fl,  we  use  dfl  =  U f=1w«.  We  take,  L 2  =  1  mm  and  L\  =  2  mm.  The 
bottom  boundary,  uq  =  {( x,y)\y  =  0,x  G  (0,Li)},  is  heated  with  heat  flux 
Sq  =  3.3  x  l0-3^from  the  initial  time,  to  =  0,  until  ts  and  insulated 
for  t  >  ts.  We  use  an  indicator  function  to  describe  the  flash- heating  of  the 
boundary, 


1  for  t  G  [t0,  ^.s] 
0  otherwise  . 


The  other  boundaries  are  insulated  throughout  the  experiment.  The  boundary 
locations  are  given  in  Figure  1.  This  experiment  can  be  described  by  the  system 
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Figure  1:  The  homogeneous  rectangle  with  boundaries. 


Cpput  —  kAu  =  0 


i = •so%,,.i(t) 

du  _ 

dp 


[u(0,£)  =  U0 


in  Q  x  (0,  T) 
on  uq  x  (0,  T) 

on  U'Li  x  (°> T) 

for  all  x  E  Q, 


(1) 


with  thermal  conductivity  k  =  3.5 x  1CP3  J/(mm  K  sec),  specific  heat  cp  =  0.75 
J/(g  K),  material  density  p  =  1.6  x  10-3  g/(mm3),  and  the  dependent  variable 
u  is  temperature  in  degrees  Kelvin.  The  values  used  throughout  this  article 
are  summarized  in  Table  1.  We  will  now  make  a  few  changes  to  the  above 
system.  In  the  partial  differential  equation,  An  may  be  written  more  generally 
as  V  •  (yf°Vu),  where  A 0  is  the  2  by  2  identity  matrix  E  in  the  present  section 
and  Section  2.  We  make  this  change  because  subsequently  we  will  replace 
the  identity  with  another  positive  definite  matrix,  which  will  be  derived  from 
homogenization  theory.  Corresponding  to  this  change,  we  must  also  change 
the  boundary  conditions,  so  that  the  boundary  conditions  still  specify  the  flux 
when  ^4°  is  not  the  identity  matrix.  Specifically,  ^  =  n  ■  Vu,  where  n  is  the 
exterior  unit  normal  vector,  is  replaced  with  =  n-  A^SJu.  For  convenience 
and  without  loss  of  generality,  we  will  also  translate  the  temperature  so  that  the 
initial  temperature,  which  we  will  take  to  also  be  the  surrounding  temperature, 
is  zero.  This  corresponds  to  the  change  u  =  u  —  uq.  We  wish  to  use  the  thermal 
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diffusivity,  a  =  ckp ,  as  a  characteristic  parameter  so  our  final  system  is 
ut  —  aV  •  (A°Vu)  =  0 

,  =  S/tw-lW 

'  “^=° 
u(x,  0)  =  0 

where  St  = 

J  CpP 

The  flash-heat  experiment  is  often  used  to  estimate  the  parameter  a.  Data 
is  typically  collected  from  the  source  boundary  (uq)  or  the  back  boundary 
(the  boundary  ui2  opposite  the  source  boundary).  In  carrying  out  parameter 
estimation,  one  minimizes  the  least  squares  distance  (or  the  Frobenius  norm) 
between  the  data  and  the  model,  where  the  Frobenius  norm  of  the  m  x  n 
matrix  A  with  entries  atJ  is 


in  Q  x  (0,  T) 
on  0J4  x  (0,  T) 

on  U'Li  x  (°’  T ) 

for  all 


(2) 


F  = 


\ 


EE' 

i=i  j= 1 


(3) 


We  will  use  the  Frobenius  norm  of  values  at  nodes  on  002  and  u> 4  and  time  points 
to  evaluate  the  difference  between  models  when  graphical  representations  are 
not  sufficient. 

For  completeness  we  will  additionally  discuss  three  modified  systems  in 
Sections  2. 1-2.3  which  are  commonly  used  to  model  heat  but  with  different 
physical  assumptions  than  those  of  (2).  In  each  of  Sections  2. 1-2.3,  we  discuss 
the  finite  element  solution  of  the  system  and  the  physical  relevance  of  the 
respective  system  in  representing  the  flash- heat  experiment.  Throughout  the 
remainder  of  this  document,  we  will  use  the  parameter  values  specified  in 
Table  1  unless  otherwise  stated. 


2.1  System  1 

The  first  modified  system  we  propose  restricts  solutions  of  (2)  to  be  in  V  = 
{u|u  €  L2(0,T;i4(D),  §  G  L2(0,T;  (iZ^fi))!)},  where  C  u;2  and  H^(Q)  = 
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Parameter 

Value 

Units 

1^0:  ts. 

[0,0.6] 

s 

T 

5 

s 

Sf 

2.75 

K  mrn/s 

a 

2.9167 

mm2/s 

Lx 

2 

mm 

L2 

1 

mm 

(U5 

{{xiD)\x  £  (1.7, 1.8), y  =  1} 

mm 

h 

l - 1 

O 

O 

l _ l 

Table  1:  Parameter  values  used  throughout  this  document,  unless  oth¬ 
erwise  mentioned.  0J5  is  an  interval  which  corresponds  to  a  small  portion 
of  the  boundary  0J2,  and  A0  is  a  2-by-2  positive  definite  matrix  arising 
from  homogenization  theory. 


{w  £  f21(ff) \w{x)  =  0,  for  x  £  CJ5}.  Here  we  denote  the  topological  dual  of  the 
space  X  as  X'.  Frequently  one  encounters  the  notation  X*  as  the  topological 
dual  of  X ,  but  we  reserve  the  use  of  *  to  represent  geometries  in  the  develop¬ 
ment  of  homogenization  results  in  Section  3.  In  the  remainder  of  this  section 
we  will  take  W5  =  {(x,y) \x  £  (1.7, 1.8),  y  =  1}.  This  can  be  incorporated 
succinctly  in  the  system 


Uf  —  aV  •  (y4°Vu)  =  0 

in  Q  x  (0,  T) 

on  W4  x  (0,  T) 

aot=° 

on  Uf=1  uji\uj5 

u  =  0 

on  W5  x  (0,  T) 

u(x,  0)  =  0. 

(4) 


We  will  use  the  finite  element  method  to  solve  (4).  The  finite  element 
method  approximates  the  infinite  dimensional  solution  of  a  partial  differential 
equation  with  a  finite  dimensional  approximation.  The  domain  (fl)  is  dis¬ 
cretized  using  the  Delaunay  triangulation.  The  finite  dimensional  solution  is 
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taken  from  the  space  of  piecewise  two  dimensional  affine  functions,  where  the 

solution  is  affine  on  each  mesh  element  (see  [15]  for  details).  Recall  the  nota- 

du 

tion  n  =  (nx,nv),  where  n  is  the  unit  outward  normal  and  — -  =  n  -  (H°V«). 

or)  ao  ' 

The  boundary  conditions  can  also  be  written  as 
du 

lw3  ~“/  lw3 


dr)Ao 

du 


du  du 

=  on  a — h  °12w- 
dx  dy 


=  0 


dyAo 

du 


du 


du 


diuo 

du 


oj  i 


OJ2 


—  —  Oil  a - «12w- 

dx  dy 


du  du 

=  021  w - 1"  °22Tp 

dx  dy 


=  0 


ui 


=  0 


a 


dr)A  o 


Cd4 


du  du 

=  a  (  021  a -  022  TT- 

dx  dy 


=  sfx[toM^- 


(5) 


id  4 


Using  Green’s  identities,  we  obtain  the  weak  form 


( p(x)utdx —  /  4>(x,  0)5/X[to  ts](f)  dx  +  a  lx7(f>-A0'Vudx  =  0 

>/  ^  J  id  4 

or 

(ut,  (j)}  -  f  (j)(x,  0)Sfl[toM  ( t)dx  +  a(V(f),  A°Vu)  =  0, 

J  id  4 

where  (4>,v)  =  j^rjrvdx.  Here  (j)  is  a  member  of  the  space  of  test  functions 
(fi)  =  {(f)  £  Hl(Q)\(j)(x)  =  0,  for  x  €  W5},  or  test  functions  which  satisfy 
the  essential  boundary  conditions.  We  approximate  u  with  un  by 


N 


uN(t,x )  =  T^\t)(j)i(x ) 


i= 1 


Substituting  this  expression  into  the  weak  form  with  test  functions  fa,  i  £ 
{1,2,...  A^}  in  the  space  of  two  dimensional  affine  functions  that  are  zero  on 
W5  (we  refer  the  reader  to  [15]  for  more  details  on  these  basis  elements  and  the 
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time  dependent  coefficients  T-  L\t)),  we  have 


d_ 

dt 


N 


E*T. 


(1M  <*■ 

i  >  9 3 


(/)j{x,0)Sfl[to,ts}(t)dx 


<i= 1 


/  CJ4 


(6) 


+  «(v^,2L0v(El=1r/1Vt))  =  o. 


We  may  factor  the  time  dependent  coefficients  T-  1  ^  (t)  of  the  basis  elements 
4>i(x)  from  the  inner  product  to  obtain 


XJ  JtTi1]  (&>  fa)  ~  SfZ[to,ta]^)  J  4>j{x,0)dx 

+  «EXi^(1)<v^,^°v^>  =  o. 


(7) 


Equation  (7)  must  be  true  for  arbitrary  j  G  {1,2,...,  N}  and  the  system  may 
be  written  as 


C 1 1  fi  (t)  +  aMjfi  (t)  =  5/%^]  (t)h ,  (8) 

where  Ci  is  an  N  x  IV  positive  definite  matrix  with  elements  =  {(j)i,(j)j), 
M\  is  an  N  x  N  positive  definite  matrix  with  elements  rri^ 1  =  (V fa,  v4°V (f>j) , 
fi  is  an  JV-vector  with  components  f-1'1  =  f  (f>i(x,0)dx  and  T\  is  an  N  col¬ 
umn  vector.  We  informally  verified  our  calculations  of  these  vectors  in  the 
case  where  A 0  =  I2  by  comparing  them  to  the  corresponding  values  used  to 
calculate  the  finite  element  method  solution  in  MatLab’s  PDE  toolbox. 

The  physical  interpretation  of  this  system  is  that  there  is  a  small  portion 
of  the  back  boundary  which  is  fixed  at  the  ambient  temperature.  This  is  not 
a  physically  feasible  situation.  In  order  to  fix  the  temperature  of  a  boundary, 
experimentalists  typically  fix  the  boundary  to  be  0  °C  by  applying  ice.  We  see 
in  Figure  2  below  that  setting  the  temperature  to  the  ambient  temperature 
does  not  give  us  the  behavior  observed  in  the  flash-heat  experiment,  so  we 
suspect  forcing  the  boundary  to  be  a  lower  temperature  (assuming  ambient 
temperature  is  above  freezing)  will  further  exacerbate  these  problems.  The 
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rapid  decrease  of  the  solution  of  (4)  to  zero  may  be  intuitively  understood  by 
imagining  being  cooled  with  ice.  The  ice  will  melt  until  the  interior  reaches 
the  temperature  of  the  ice,  and  the  resulting  water  (from  the  melted  ice)  is 
related  to  the  loss  of  energy  over  time.  We  could  make  005  smaller,  but  the 
number  of  mesh  points  required  would  increase  significantly  to  do  this  and 
the  computational  costs  could  become  exorbitant.  So  to  use  this  model  we 
would  have  to  alter  the  experimental  set-up  to  accommodate  setting  a  small 
part  of  the  boundary  to  a  fixed  temperature.  The  inefficacy  of  this  system  to 
capture  the  behavior  of  the  flash-heat  experiment  is  demonstrated  in  Figure  2 
where  we  compare  the  solution  of  (4)  to  the  solution  of  (2)  with  A 0  =  E-  In 
Figure  2(a)  and  (b),  we  see  that  the  temperature  of  the  midpoint  of  both  the 
front  and  back  boundaries  go  to  zero  in  time.  It  is  also  important  to  note  that 
(4)  does  not  capture  the  shape  of  temperature  on  the  back  boundary  U2  at 
t  =  1  as  we  see  in  Figure  2(c).  We  chose  time  t  =  1  for  Figure  2  because  it 
demonstrates  the  lack  of  smoothness  in  the  solution  though  it  is  still  relatively 
close  the  the  solution  of  the  flash-heat  experiment  as  compared  to  the  solution 
at  t  =  5  (see  Figure  2(a)).  Because  our  colleagues  at  NASA  Langley  Research 
Center  do  not  currently  have  experiments  which  fix  the  temperature  of  the 
boundary  of  a  specimen,  we  will  not  use  (4)  in  the  remainder  of  this  work. 


Temperature 
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Temperature  at  *  =  1,  y  =  1  Temperature  at  *  =  1,  y  =  0 


Back  boundary  at  t  =  1  sec 


Figure  2:  (a)  The  solutions  of  (2)  and  (4)  labeled  Flash  Experiment  and 
System  1,  respectively,  at  the  midpoint  of  the  back  boundary;  (b)  The 
solutions  of  the  (2)  and  (4)  at  the  midpoint  of  the  source  boundary;  (c) 
The  solutions  of  the  (2)  and  (4)  on  the  back  boundary  at  t  —  1. 
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2.2  System  2 

The  second  system  generalizes  the  operator  —V  •  A°V  with  the  operator  —V  • 
A°V  +  A  for  A  >  0.  The  difference  between  the  solution  of  this  system  and 
the  solution  of  (2)  will  depend  on  the  magnitude  of  A  and  clearly  there  is  no 
difference  for  A  =  0.  We  will  use  A  =  10-2  to  study  the  behavior  of  this  system 
and  will  further  discuss  the  choice  of  A  at  the  conclusion  of  this  section.  The 
modified  system 


ut  —  aV  •  (A°Vu)  +  Xu  =  0 

in  17  x  (0,  T) 

“e£.  =  M 

on  u: 4  x  (0,  T) 

Co  =  0 

on  Uf=1  cut  x  ( 

(9) 


u(x,  0)  =  0, 


has  a  Xu  term  in  the  partial  differential  equation.  The  Xu  term  represents 
loss  on  the  entire  specimen  in  the  direction  orthogonal  to  the  two  dimensional 
specimen,  while  there  is  no  loss  on  the  boundary  of  the  specimen. 

First,  we  will  consider  the  numerical  solution  of  (9).  Again  using  the  finite 
element  method  and  the  steps  detailed  in  Section  2.1,  we  begin  by  finding  the 
weak  form  of  (9)  given  by 

[ip  ( ut  +  Xu)  +  aSJip  ■  A0  Vu]  dx  —  Sfl[t(hts ]  /  ip(x,  0)  dx  =  0, 

y  ^  J  W4 

with  test  functions  ip  £  i/1(fl).  We  chose  test  functions  from  Hl(Q.)  because 
there  are  no  essential  boundary  conditions  in  this  formulation.  We  replace  u 
with  the  finite  dimensional  approximation  ul  defined  by 

L 

uL(x,y,t )  =  ] 'TT-2\t)ipi(x,y ), 

i= 1 

where  T”  (t)  are  the  time-dependent  coefficients  of  the  two  dimensional  piece- 
wise  affine  basis  elements  ipi  which  are  described  in  more  detail  in  [15].  Using 
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similar  calculations  as  for  (4),  we  have 


(10) 


=  Sfx[t0lts]  /  i/>j(x,0)dx, 


where  (if),  v )  =  il)v  dx.  We  write  this  system  as  a  system  of  ordinary  differ¬ 
ential  equations 


c2£  T2(t )  +  aM2T2(t)  =  Sfl[toM(t)f2,  (11) 

where  C2  is  an  L  x  L  positive  definite  matrix  with  entries  c\  ■  =  (?/>j,  M2  is 

(2)  a  J 

anlxL  positive  definite  matrix  with  entries  m\-  =  ^ (^ ,  ^-}  +  (V'i/’i,  AoVV’j) , 
f2  is  an  L- vector  with  components  =  f  V’ i(x >  0)dx  and  T2  is  an  L- vector. 

As  we  can  see  in  Figure  3  below,  the  solution  of  (9)  is  more  similar  to  the 
solution  of  (2)  than  is  the  solution  of  (4).  There  will  certainly  be  a  decrease 
in  the  difference  between  the  solutions  of  (2)  and  (9)  toward  zero  as  A  goes 
to  zero.  Like  the  solution  to  (4),  solutions  of  this  limit  system  will  approach 
zero,  though  it  will  take  much  longer  than  the  solution  of  (4).  The  decrease  of 
solutions  to  (9)  toward  zero  is  most  clear  in  Figure  3(a)  and  (b).  An  advantage 
of  this  system  is  that  the  shape  of  the  temperature  at  the  boundaries  U14  and 
uj2  is  similar  in  shape  to  the  solution  of  (2).  For  instance,  in  Figure  3(c)  where 
t  =  5  which  is  the  time  that  gives  the  largest  difference  between  the  two  models 
on  the  time  interval  [0,5],  the  temperature  on  the  back  boundary  is  constant 
across  the  back  boundary  for  both  the  solution  of  (2)  and  (9).  If  we  take  A 
small  and  the  final  time  T  not  too  large,  (9)  could  be  a  good  approximation 
of  (2).  For  A  =  10~5,  the  Frobenius  norm  of  the  difference  between  uflash  the 
solution  of  (2)  and  the  solution  of  (9)  U  at  the  back  boundary  u2  and  source 
boundary  UJ4  are  nearly  equal  and  very  small.  For  instance,  the  Frobenius 
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norm  taken  on  uq  and  t  £  [0,  5]  is 
41  101 

EE  (Ufa,  0,  tj)  -  u&ashfa,  0,  tj)f  =  0.0028, 

i=  1  j=  1 

where  xt  =  and  tj  =  ■  These  advantages  in  the  example  where  the 

partial  differential  equation  is  defined  on  a  homogeneous  rectangle  (Q)  with 
isotropic  flow  ( A 0  =  I2)  lead  us  to  choose  (9)  with  A  =  10-5  to  model  the 
flash-heat  experiment  in  Sections  3  and  4  where  we  relax  the  assumptions  on 
the  domain  of  the  partial  differential  equation  and  consider  anisotropic  flow. 


Temperature 
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Temperature  at  x  =  1,  y  =  1  Temperature  at  X  =  1,  y  =  0 


(a) 


(b) 
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0.2 
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Back  boundary  at  t  =  5  sec 
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- Flash  Experiment 

-i —  System  2 


0 


0.5  1 

Position  on  W2 

(c) 


1.5 


Figure  3:  (a)  The  solutions  of  (2)  and  (9),  labeled  Flash  Experiment 
and  System  2,  respectively,  with  A  =  0.01  at  the  midpoint  of  the  back 
boundary;  (b)  The  solutions  of  the  (2)  and  (9)  with  A  =  0.01  at  the 
midpoint  of  the  source  boundary;  (c)  The  solutions  of  the  (2)  and  (4) 
with  A  =  0.01  on  the  back  boundary  at  t  —  5. 
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2.3  System  3 

The  most  physically  reasonable  mathematical  model  of  the  flash-heat  experi¬ 
ment  replaces  the  zero-flux  boundary  conditions  of  (2)  with  Robin  boundary 
conditions.  Robin  boundary  conditions  relax  the  assumption  of  perfectly  in¬ 
sulated  boundaries  by  assuming  Newton  cooling  occurs  on  the  boundaries.  In 
this  section  we  consider  two  cases,  the  first  has  Robin  boundary  conditions 
only  on  co^  C  0J2 


ut  —  aV  •  (yl0Vit)  =  0 

in  12  x  (0,  T) 

“4:. = 

on  CJ4  x  (0,  T) 

Co=° 

on  Uf=1  uji\uj5  x  (0,  T) 

<UA0= 

on  W5  x  (0,  T) 

w(x,  0)  =  0, 


(12) 


where  W5  =  {(x,y)\x  G  (1.7, 1.8), y  =  1}.  In  the  second  case  we  include  Robin 
boundary  conditions  in  the  entire  boundary  by  incorporating  a  —Xu  term  in 
all  of  the  boundary  conditions 


ut  —  «V  •  (R°Vn)  =0  in  11  x  (0,  T) 

<  =  SfX[toM  ( t )  -  Xu  on  cu4  x  (0,  T) 

=  ~^u  on  u*3=i UJi  X  ('0,  T ) 

u(x,  0)  =  0. 

It  is  important  to  recall  from  earlier  in  Section  2  that  u  =  u  —  uq  where  u  is 
the  temperature  and  uq  is  both  the  initial  temperature  and  the  temperature 
surrounding  the  specimen.  To  solve  (12)  and  (13)  numerically  with  the  finite 
element  method,  we  begin  with  the  weak  form  of  (12)  and  (13) 


f  (~fut  +  aVy  •  A°X7u)  dx  -  Sfl[tojs](t)  f  j(x,0 )dx 

w  J  (J4 


(14) 


+  X  'yu  ds  =  0 

Jdi 1 
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with  7  £  i71(J7),  where  <911  is  u> 5  for  (12)  and  d£l  for  (13).  We  approximate  u 
with  up  defined  by 

p 

uP{x,y,t)  =  EE^biO,!/), 

i=  1 

where  7 \  are  two  dimensional  piecewise  affine  basis  elements  with  time  depen- 

/o\ 

dent  coefficients  T-  (t)  described  in  more  detail  in  [15].  After  calculations 
similar  to  those  reported  in  Section  2.1,  we  obtain  the  expression 


%—  1 


2—1 


(15) 


P 

sfT[to,ts] (*)  [  7 3 (*>  0)  dx  +  A  f  E  Ti3) (thjH  ds, 

•J  Qi 4  j 


for  arbitrary  j  £  {1,  2  . . .  P}.  The  resultant  system  can  be  written  as 

c3jt  ut)  +  W  +  AH)  f3(f)  =  5/X[t0)t.](t)/3,.  (16) 

(o\  /o\ 

where  C3,M3  and  D  are  P  x  P  matrices  with  entries  c-  -  =  (7*,  7j).  = 

(7j,A°7j),  and  dij  =  ds,  respectively,  and  /3  is  a  P- vector  with  com¬ 

ponents  /j3)  =  /W4  7j(x)  da- 

The  boundary  conditions  in  (2)  correspond  to  zero  flux  on  the  boundaries, 
which  physically  represents  perfectly  and  flawlessly  insulated  boundaries.  The 
system  (12)  relaxes  the  flawless  assumption.  The  Robin  boundary  conditions 
imposed  on  W5  in  (12)  would  correspond  to  a  flaw  in  the  insulation  at  W5  which 
leaches  heat  with  convection  heat  transfer  coefficient  A.  The  system  (13)  also 
relaxes  the  assumption  that  the  insulation  is  perfect  with  the  additional  Xu 
term  in  all  of  the  boundary  conditions  corresponding  to  an  insulating  mate¬ 
rial  which  has  convection  heat  transfer  coefficient  A.  The  somewhat  natural 
physical  interpretation  of  these  two  systems  is  an  advantage  as  compared  to 
the  modified  systems  we  considered  in  Sections  2.1  and  2.2. 

In  Figure  4  we  see  that  with  A  =  0.1  the  solution  of  (12)  over  time  on  the 
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source  and  back  boundary  is  quite  close  the  solution  of  (2).  The  solution  of 
(13)  with  A  =  0.1  and  A0  =  E  decreases  toward  zero  much  more  than  the 
solution  of  (13)  as  is  clearly  demonstrated  at  the  midpoints  of  0J2  and  uq  in 
Figures  5(a)  and  (b),  respectively.  However,  if  we  decrease  A  by  a  few  orders 
of  magnitude,  the  solution  of  (13)  should  possess  a  decrease  similar  to  the 
decrease  seen  in  the  solution  of  (12)  depicted  in  Figures  4(a)  and  (b). 

We  hope  to  use  the  models  developed  in  this  effort  in  the  future  to  model 
experiments  which  will  be  carried  out  at  NASA  Langley  Research  Center.  It 
is  a  reasonable  assumption  that  the  insulating  material  used  in  performing  the 
flash-heat  experiment  will  not  be  damaged-the  assumption  corresponding  to 
(12),  especially  as  we  hope  this  work  may  be  used  to  characterize  damage  in 
porous  materials.  However,  the  assumption  that  the  insulation  used  in  practice 
is  not  perfect  and  allows  some  exchange  of  heat  between  the  specimen  and  the 
surroundings  is  not  only  feasible  but  inevitable.  The  Frobenius  norm  of  the 
difference  between  uflash,  the  solution  of  (2),  and  U,  the  solution  of  (13),  on 
the  source  boundary  uq  is  sufficiently  small  for  A  =  10~'5.  More  precisely, 
|| U  —  uflash||F  on  UJ4  and  t  £  [0, 5]  is  given  by 

=  0.0085, 

where  Xj  =  ^  and  tj  =  .  The  physical  relevance  of  (13)  motivates  us 

to  use  (13)  with  A  =  10-5  in  addition  to  (9)  with  A  =  10-5  to  model  the 
flash-heat  experiment  in  the  subsequent  sections. 
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Temperature  at  x  =  1,  y  =  1 


(a) 


Temperature  at  x  =  1,  y  =  0 


(b) 
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Back  boundary  at  t  =  5  sec 


H — I — I — I — I— I — I — I— I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — | — I — I — | — | — | — I- 


0 


- Flash  Experiment 

-t —  System  3 


0.5  1  1.5 
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Figure  4:  For  dQ.  =  cu5  and  A  =  0.1  (a)  The  solutions  of  (2),  labeled 
Flash  Experiment,  and  (12),  labeled  System  3,  at  the  midpoint  of  the 
back  boundary;  (b)  The  solutions  of  the  (2)  and  (12)  at  the  midpoint  of 
the  source  boundary;  (c)  The  solutions  of  the  (2)  and  (12)  on  the  back 
boundary  at  t  =  5. 
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Temperature  at  x  =  1,  y  =  1  Temperature  at  x  =  1,  y  =  0 


(a) 


(b) 


Back  boundary  at  t  =  5  sec 
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Figure  5:  For  dil  =  dil  and  A  =  0.1  (a)  The  solutions  of  (2),  labeled 
Flash  Experiment,  and  (13),  labeled  System  3,  at  the  midpoint  of  the 
back  boundary;  (b)  The  solutions  of  the  (2)  and  (13)  at  the  midpoint  of 
the  source  boundary;  (c)  The  solutions  of  the  (2)  and  (13)  on  the  back 
boundary  at  t  —  5. 
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3  Homogenization 

Now  that  we  have  chosen  two  systems,  (9)  with  A  =  10-5  and  (13)  with 
A  =  1CH5,  which  model  the  flash-heat  experiment  reasonably  well  on  a  ho¬ 
mogeneous  domain,  we  will  use  these  two  systems  to  model  the  flash-heat 
experiment  on  more  complicated  domains  than  the  homogeneous  rectangle  Q. 
We  developed  a  method  to  model  the  flash-heat  experiment  on  porous  do¬ 
mains  in  [6].  Though  these  simulations  were  useful  in  detecting  damage,  they 
were  too  computationally  intensive  for  the  more  sophisticated  parameter  esti¬ 
mation  procedures  needed  to  characterize  damage.  Here  we  will  use  (9)  and 
(13)  along  with  random  geometries  to  model  the  flash- heat  experiment  in  a 
porous  domain  and  compare  them  to  limit  partial  differential  equations  which 
are  derived  from  homogenization  theory.  In  this  formulation  A0  is  no  longer 
1-2  and  the  random  complicated  geometry  is  replaced  with  a  less  complicated 
domain  with  anisotropic  flow.  We  will  focus  on  values  obtained  on  the  source 
boundary  ui 4  when  we  compare  the  subsequent  models. 

We  will  consider  a  random  geometry  Q  (depicted  in  Figure  6),  which  is 
composed  of  H  \  j  ,  where  are  randomly  placed  pores,  <912*  is  the 
boundary  of  the  ith  pore  and  c u*  are  the  same  as  in  the  previous  sections.  For 
each  of  the  systems  (9)  and  (13)  we  chose  on  the  homogeneous  rectangle  fi,  we 
must  now  pose  the  problem  on  H.  We  will  call  the  dependent  variable  of  this 
system  urand,  where  ‘rand’  refers  to  the  random  domain.  Here  it  is  worthwhile 
to  note  that  we  have  performed  the  transformations  detailed  in  the  beginning 
of  Section  2  to  all  of  the  partial  differential  equations  in  this  section. 

The  system  corresponding  to  (9)  on  the  random  domain  fl  is  given  by 


u$and  -  aAurand  +  Aurand  =  0 

on 

f l  X 

h 

0 

a  J^rand  = 

0 

on 

1  1  Tlr 

idfiiX  (0,T) 

a  J^rand  = 

0 

on 

uf= 

1  Ui  x  (0,  T) 

a  J^rand  = 

-S' 

0 

H 

on 

o’ 

X 

-+ 

3 

itrand(x,  0)  = 

=  0, 

< 


(17) 
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co2 


Figure  6:  A  typical  porous  domain  with  random  pores  (enlarged  view) 


and  the  system  corresponding  to  (13)  on  the  random  domain  ft  is  given  by 
'  u\and  -  aAurand  =  0  on!lx  (0,  T) 

a|_nrand  _  q  Gn  QQ.  x  (0>  T) 

<  a^uiand  =  —Xuiand  on  uf=1  coi  x  (0 ,T) 

a|_nrand  =  Sfl[0)ts]  (t)  -  Aurand  on  cu4  x  (0,  T) 
urand(x,  0)  =  0. 

We  will  use  these  two  systems  on  random  geometries  to  simulate  data 
that  one  might  expect  from  the  flash-heat  experiment  performed  on  porous 
specimens.  We  generate  these  geometries  and  solve  the  partial  differential 
equations  using  the  methods  developed  in  [6]. 

In  order  to  apply  homogenization  theory  to  (17)  and  (18),  we  will  use  ge¬ 
ometries  with  enough  pores  to  suppose  that  ft  has  a  periodic  structure  (though 
it  may  be  that  the  physical  specimens  are  better  modeled  with  periodic  pores 
as  the  pores  in  composite  materials  are  often  the  result  of  sinusoidal  manu¬ 
facturing  processes).  In  other  words,  we  suppose  that  the  pores  (or  holes)  are 
periodically  distributed  with  a  period  e,  where  £  is  a  small  parameter  that  we 
let  go  to  zero  in  the  limit.  This  is  the  framework  of  the  periodic  homogeniza¬ 
tion  theory  which  is  explained  in  more  detail  in  [11].  To  do  so,  we  introduce 


(18) 
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a  reference  cell  (or  domain)  Y .  For  our  purpose,  let  us  take  as  Y  the  origi¬ 
nal  (homogeneous)  rectangle  Q.  Let  B  =  U^B*  be  a  set  of  N  open  subsets 
strictly  included  in  Y  such  that  Bj  n  Bj  =  0  for  i  %  j. 

Denote  by  r(eB)  the  set  of  all  translated  images  of  eB  of  the  form  e(k£+B), 
k  €  Z2,  kI  =  (ki£i,  K2O),  so  the  set  r(eB)  represents  the  periodic  pores  in 
R2.  Let  B£  be  the  set  of  the  holes  contained  in  fL  With  the  above  choice  of 


Y,  and  taking  for  instance  e  = 


r/¥j 


with  n£  — >  +00,  no  hole  from  Be  will 


intersect  the  boundary  dQ.  Here  [~-J  is  the  nearest  integer  function.  We  will 
set 

n£  =  n\B£. 

By  this  construction,  the  physical  domain  is  periodically  perforated  with 
holes  of  size  of  the  same  order  as  the  period.  We  are  essentially  approximating 
the  random  geometry  Q  with  the  periodic  geometry  fl£. 

We  will  use  the  following  notation: 

•  Y*  =Y\B, 

■  X\  the  topological  dual  space  of  X 


e  = 


\Y* 

IFi 


-,  the  proportion  of  the  material  in  the  cell  Y, 


■  |w|  =  the  Lebesgue  measure  of  any  open  set  a ;, 

•  — r  [  <p(x)  dx,  the  mean  value  of  <p  on  the  set  u. 

Iwl  JiO 

Observe  that  by  construction,  0  is  also  the  proportion  of  the  material  in 
for  any  e  >  0  and  the  percent  porosity  =  (1  —  0)  x  100%.  We  are  now 
prepared  to  consider  our  two  systems  on  He.  Using  fle  to  approximate  H  in 
(17)  we  have 


uf  —  aAu£  +  A u£  =  0  in  Qe  x  (0,  T) 

on  dB£  U?=1  ojt  x  (0,  T ) 
on  W4  x  (0,  T) 

u£(x,  0)  =  0, 


a^u£  =  0 

=  5/X[o,ts](i) 


(19) 
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and  in  (18) 

uf  —  aAu£  =  0  in  x  (0,  T) 

=  0  on  dB£  x  (0 ,T) 

<  aJ^u£  =  -A u£  on  uf=1  ui  x  (0,  T)  (20) 

ot-^u£  =  Sfl[ 0:ts]  ( t )  -  A u£  on  u4  x  (0,  T) 

ne(.x,  0)  =  0. 

Both  (19)  and  (20)  have  unique  solutions  in  the  Banach  space  W£  =  {w|u  G 
L2(0,  T;  _ff1(He)),  G  L2(0,  T;  (_ff1(fl£))t)}.  One  is  then  allowed  to  pass  to 
the  limit  in  (19)  and  (20)  to  obtain  the  limit  homogenized  systems  (for  details, 
we  refer  the  reader  to  [10,  11,  12]).  In  particular,  these  references  contain 
proofs  that 

u£  — ^  0U  weakly  in  L2(0,T;  .ft1  (H)),  (21) 

where  U  is  characterized  by  the  unique  solution  of  homogenized  problems  and 
u£  is  the  zero  extension  of  u£  from  to  the  whole  domain  A. 

The  limit  system  corresponding  to  (19)  is  given  by 

'  9Ut  -  aV  ■  (A°U)  +  X9U  =  0  in  x  (0 ,T) 

a  Au  =  0  on  U?i  uji  x  (0,  T ) 

°vAo  *-!  v  ’  7  (22) 

=  SfX\ o.t.] (*)  on  W4  x  (°> T) 

U(x,  0)  =  0, 

while  the  limit  system  corresponding  to  (20)  is  given  by 

eut  -  aV  ■  (A°U)  =  0  in  11  x  (0,  T) 

aJ^-  =  -\U  on  uL  w,  x  (0,T) 

°vAo  t-i  V  >  (23) 

=  Sfl[0M(t )  -  At/  on  x  (0,T) 
k^(*,0)  =  0. 
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In  both  limit  systems  the  homogenized  operator  A 0  has  constant  coefficients 
and  is  expressed  in  terms  of  the  following  cell  problems  defining  the  “corrector” 


3  HOMOGENIZATION 
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functions  xi  and  X2- 


—A Xj  =  0  for  j  =  1,2  in  Y* 
Xj  is  Y  periodic 
0 

(Xj  -Vj)  =  0  on  dtti 

M*y(x)  =  0. 

Then  the  homogenized  matrix  ^4°  =  (oF)  is  defined  by 


o  a  1  I  dxi a  o  1 
an  =  ^-jT7i  -z-dy,  a12  =  -  — 


1 


\Y 


lY* 


dy  i 


<9x2 

\Y\  Jy*  d'!J\ 
9X2 


0  1  /'  <9X1  ,  o  _n  1 

“21  “  “ W\  Jy  w/y'  a^-e^W\  Jy  % 2 


dy, 

dy. 


(24) 


(25) 


It  is  important  to  note  that  <9fh  in  (24)  refers  to  the  boundaries  of  the  pores 
of  the  reference  cells,  depicted  in  Figure  7. 


Ib.,  s=l:2 
'  2 


£L.  s=1M 

2 


Reference  Cel! 

Figure  7:  The  reference  domain  and  the  limit  domain  with  £  = 


\  and 
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It  can  be  shown  that  the  error  estimate  (distance  between  u£  and  6U  in 
the  L2(0,  T;  H1  (Ll))-noicm)  is  of  order  of  yfi,  which  completely  justifies  the  ho¬ 
mogenization  procedure  if  e  is  sufficiently  small,  or  equivalently  if  the  number 
n£  of  holes  is  sufficiently  large  since  n£  ~  — r.  Using  the  linear  trace  operator, 
7  :  Hl(U)  — y  L2(ujR)  and  arguments  similar  to  those  in  [8],  we  get  the  result 
that  7 (u£)  converges  weakly  in  L2(ujR)  to  7 (U). 

Recalling  the  convergence  of  u£  to  U,  we  created  numerical  simulations  to 
compare  U  to  u£  and  ttrand  on  0J4 .  We  used  methods  developed  in  [6]  and  [16]  to 
generate  the  random  geometries  P  and  Y* .  We  used  MatLab’s  PDE  toolbox, 
as  in  [6],  to  then  solve  the  partial  differential  equations  on  these  domains.  To 
calculate  U,  we  used  the  finite  element  schemes  detailed  in  Sections  2.2  and 
2.3. 

For  each  simulation,  letting  N  be  the  number  of  pores  in  the  reference  cell, 
and  nr  the  number  of  pores  in  the  random  geometry,  we  might  suppose  that 
the  e  corresponding  to  the  f le  which  approximates  Q  is  given  by 


£  = 


(26) 


with  [~-J  representing  the  nearest  integer  function.  As  we  see  in  (26),  e  decreases 
as  N  decreases.  This  leads  to  a  subtlety  in  choosing  N .  We  would  like  Y* 
to  capture  the  random  nature  of  P  while  still  containing  a  sufficiently  small 
number  of  pores  N  to  ensure  £  is  small.  In  the  simulations  presented  here,  we 
take  N  =  2  and  use  (26)  above  to  calculate  u£ . 


4  Results 

4.1  Results  using  systems  corresponding  to  (9) 

In  the  previous  section  we  developed  (17)  and  (with  solution  urand)  to  model 
the  flash- heat  experiment  in  porous  material  and  introcuded  (19)  (with  solu¬ 
tion  u£ )  to  approximate  (17).  Finally,  we  introduced  homogenization  theory  to 
approximate  (19)  with  (22)  (with  solution  U).  In  this  section  we  will  present 
graphical  results  corresponding  to  the  three  systems  and  summarize  the  re- 
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suits  of  200  simulations  of  the  three  systems.  Here,  we  will  use  the  parameter 
values  listed  in  Table  1  with  A  =  10-5.  In  the  next  section,  we  will  present 
similar  results  but  for  the  more  physically  relevant  system  corresponding  to  the 
flash-heat  experiment  on  the  perforated  domain  represented  by  system  (18). 

We  performed  simulations  of  three  different  porosity  levels,  10%  (depicted 
in  Figure  8),  5%  (depicted  in  Figure  9),  and  2%  (depicted  in  Figure  10).  The 
random  geometries  used  for  the  examples  are  featured  in  Figures  8— 10(a), 
which  corresponds  to  used  in  (17)  to  solve  for  urand.  The  approximations  of 
H,  He,  are  presented  for  each  level  of  porosity  in  Figures  8-10(b).  Figures  8- 
10(c)  give  the  solution  of  urand,  u£  and  U  at  y  =  0,  t  =  0.1  on  CJ4.  We  note 
that  U  seems  to  capture  the  ‘average’  behavior  of  urand,  and  u£.  There  are 
several  things  that  we  note  in  Figures  8-10.  For  instance  the  solution  U  to  (22) 
underestimates  -urand  at  the  x  location  where  there  are  pores  near  004  in  the 
random  geometry  f2.  For  instance,  for  8  =  0.9  in  Figure  8(a)  we  see  large  pores 
near  W4  at  x  =  0.4  and  x  =  1.6  which  correspond  to  two  peaks  in  itrand  where 
itrand  deviates  above  U  at  x  =  0.4  and  x  =  1.6  in  Figure  8(c).  However,  in 
Figure  8(b)  there  are  no  pores  in  near  U4  and  U  appears  to  overestimate  u£ 
along  the  entire  interval  x  6  [0,  2]  in  Figure  8(c).  In  Figure  9(a)  with  8  =  0.95 
there  are  no  large  pores  near  UJ4  in  fl  and  accordingly  there  are  no  peaks  in 
itrand  away  from  U  though  there  are  a  few  small  bumps  of  «rand  above  U  at 
x  =  0.3  and  between  x  =  0.75  and  x  =  1  corresponding  to  a  small  pore  near  U4 
and  a  cluster  of  large  pores  about  0.05  mm  away  from  U4,  respectively.  Again 
there  are  no  pores  in  U£  near  U4  in  Figure  9(b)  and  u£  is  overestimated  by  U 
in  Figure  9(c).  In  Figure  10(b),  where  8  =  0.98,  there  are  pores  very  close  to 
u>4  in  and  unlike  in  Figures  8(c)  and  9(c)  U  underestimates  u£  where  the 
pores  are  close  to  U4  in  IT.  In  Figure  10(a)  there  is  cluster  of  pores  near  U4  at 
x  =  1.6,  which  corresponds  to  a  spike  in  «rand  above  U  in  Figure  10(c).  It  is 
also  worthwhile  to  note  that  as  8  increases  from  0.9  to  0.98,  ttrand  approximates 
better  u£  which  approximates  better  U  in  Figures  8-10(c).  This  suggests  that 
at  lower  levels  of  porosity  homogenization  gives  us  a  better  approximation  of 
urand  and  u£  than  at  higher  levels  of  porosity. 
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Random  Sample,  6  —  0.9 
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Figure  8:  An  example  simulation  at  9  =  0.9,  or  equivalently  10%  poros¬ 
ity  (a)  The  random  geometry  hi;  (b)  The  geometric  approximation  of  11 
assumed  in  homogenization  theory,  fl£  where  e  =  0.14;  (c)  ■urand(;r,  0,  t), 
u£(x,0jt),  and  U(x,Q,t )  for  t  =  0.1  sec  on  (u4.  For  a  better  repre¬ 
sentation  of  the  behavior  of  the  three  systems  over  time  on  the  source 
boundary  (oq),  see  the  corresponding  movie  by  clicking  here. 
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Random  sample,  9  =  0.95 
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Figure  9:  An  example  simulation  at  9  =  0.95,  or  equivalently  5%  poros¬ 
ity  (a)  The  random  geometry  hi;  (b)  The  geometric  approximation  of 
assumed  in  homogenization  theory,  where  e  =  0.17;  (c)  urand(x,  0,  t), 
u£(x,0jt),  and  U(x,0,t)  for  t  =  0.1  sec  on  ui4.  For  a  better  repre¬ 
sentation  of  the  behavior  of  the  three  systems  over  time  on  the  source 
boundary  (024),  see  the  corresponding  movie  by  clicking  here.. 
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Random  Sample,  9  =  0.98 
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Figure  10:  An  example  simulation  at  6  =  0.98,  or  equivalently  2%  poros¬ 
ity  (a)  The  random  geometry  £4;  (b)  The  geometric  approximation  of  £4 
assumed  in  homogenization  theory,  £4e  where  e  =  0.17;  (c)  urand(;c,  0,  t), 
ue(x,0,t),  and  U(x,0,t)  for  t  =  0.1  sec  on  oq.  For  a  better  repre¬ 
sentation  of  the  behavior  of  the  three  systems  over  time  on  the  source 
boundary  (uq),  see  the  corresponding  movie  by  clicking  here.. 
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Model 

e  =  0.9 

9  =  0.95 

9  =  0.98 

Mean  time 

^rand 

47  sec 

51  sec 

68  sec 

u£ 

23  sec 

28  sec 

42  sec 

U 

4.0  sec 

3.8  sec 

4.1  sec 

Mean  mesh  size 

^rand 

7855 

8571 

10111 

u£ 

4882 

5780 

7484 

U 

355 

355 

355 

Table  2:  In  the  first  three  rows,  the  last  three  columns  contain  the 
average  time  for  the  200  simulations  to  calculate  the  respective  solution. 
The  last  three  rows  contain  the  average  mesh  size  over  200  simulations. 


We  carried  out  200  simulations  at  each  porosity  level  ( 9  =  0.9, 0.95  and 
0.98)  similar  to  those  used  to  generate  Figures  8-10.  For  each  simulation, 
we  selected  a  random  domain  Q  and  a  random  reference  cell  to  create  the 
periodically  perforated  domain  fl£.  We  used  these  domains  to  solve  for  urand, 
u£  and  U,  as  we  discussed  in  Section  3.  We  recorded  the  total  time  to  solve 
each  system,  the  size  of  the  mesh  used  in  the  associated  finite  element  method 
and  the  Frobenius  norm  of  the  difference  between  solutions  on  CJ4.  In  Table  2 
we  report  the  average  time  to  solve  each  of  the  systems  and  the  average  mesh 
size  used  to  discretize  the  respective  domain.  It  is  clear  that  of  the  three 
systems  (17)  is  the  most  computationally  intensive  with  a  significantly  larger 
mesh  than  the  mesh  used  to  solve  for  u£.  The  limit  system  (22)  is  the  least 
computationally  intensive;  its  mesh  is  an  order  of  magnitude  smaller  than  the 
meshes  for  the  other  two  systems.  Also  the  computing  time  required  to  solve 
(22)  for  U  is  five  to  ten  times  less  than  the  time  to  solve  (19)  for  u£  and  ten  to 
fifteen  times  less  time  to  solve  (17)  for  urand.  This  supports  the  use  of  the  limit 
system  (22)  as  a  more  computationally  suitable  model  than  (17)  and  (19). 

The  limit  system  (22)  seems  to  be  a  good  approximation  of  (17)  and  (19). 
We  evaluate  this  statement  by  considering  the  Frobenius  norm.  Here  the 
Frobenius  norm  is  taken  at  time  points  tj  =  for  j  E  {1,2,  ...101}  and 
nodes  x%  =  ^  for  i  E  {1,2, . .  .41}  along  the  source  boundary  uq.  For  each 
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9  =  0.9 

9  =  0.95 

9  =  0.98 

Mean  |  urand  —  U\\p 

1.21 

0.630 

0.268 

St.  Dev.  ||wrand  -  U\\F 

0.246 

0.106 

0.046 

Mean  \\u£  —  U\\F 

1.363 

0.669 

0.265 

St.  Dev.  | \u£  —  U\\f 

0.250 

0.128 

0.043 

Table  3:  The  first  and  third  row  contain  the  mean  Frobenius  norm  of 
the  difference  of  the  solutions  for  the  200  simiulations.  The  second  and 
fourth  row  contain  the  standard  deviation  of  the  Frobenius  norm  of  the 
difference  of  the  solutions  for  the  200  simulations. 


simulation,  we  calculated  the  Frobenius  norm  of  \\U  —  «rand||^ 

II U  -  urand||F  =  y/gii  £5=i  {U{xi,  0,  tj)  -  u^(Xi,  0,  tj))\  (27) 

and  \\U  —  u£\\f-  We  report  the  average  and  standard  deviation  of  these  quan¬ 
tities  for  each  porosity  level  in  Table  3.  As  6  increases,  U  better  approximates 
both  itrand  and  u£.  However,  it  is  not  clear  whether  U  better  approximates  u£ 
versus  ttrand. 
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4.2  Results  using  systems  corresponding  to  (13) 

We  next  report  on  computation  similar  to  those  in  the  previous  section,  but 
rather  than  using  (17),  (19)  and  (22)  as  models,  we  will  use  (18),  (20)  and  (23). 
In  each  figure  urand  is  the  solution  of  (18),  u£  is  the  solution  of  (20)  and  U  is 
the  solution  of  (23).  We  will  consider  the  three  different  porosity  levels,  10% 
(or  9  =  0.9,  depicted  in  Figure  11),  5%  (or  9  =  0.95,  depicted  in  Figure  12), 
and  2%  (or  9  =  0.98,  depicted  in  Figure  13).  The  results  presented  here  are 
for  A  =  10~5  though  we  carried  out  calculations  for  A  =  0  (the  perfectly 
insulated  model)  and  obtained  similar  results.  The  random  geometries  used 
for  the  examples  is  featured  in  Figures  ll-13(a),  which  corresponds  to  P  used 
in  (18).  The  approximations  of  P,  £ls,  are  presented  for  each  level  of  porosity 
in  Figures  11-1 3(b).  Figures  11-1 3(c)  give  the  solutions  urand,  u£  and  U  at 
y  =  0,  t  =  0.1  over  all  values  of  x  on  U4  or  the  source  boundary.  In  this  case, 
we  see  that  U  seems  to  capture  the  ‘average’  behavior  of  urand  and  that  of  u£ . 
We  observe  similar  trends  as  discussed  in  the  previous  section,  specifically  that 
U  overestimates  when  there  are  large  pores  near  oj 4  and  underestimates  when 
there  are  not  large  pores  near  0J4.  In  Figures  11  with  9  =  0.9  we  see  that  the 
cluster  of  pores  near  UJ4  in  Figure  11(a)  corresponds  to  U  underestimating  •urand 
for  x  greater  than  1.25  in  Figure  11(c).  There  are  pores  near  0J4  in  Figure  11(b) 
which  corresponds  to  periodic  peaks  of  u£  above  U  in  Figure  11(c).  There  are 
pores  near  0J4  in  in  Figure  12(b)  but  they  are  very  small  and  U  overestimates 
u£  in  Figure  12(c).  As  should  be  expected  for  these  systems,  the  difference 
between  the  three  solutions  seems  to  decrease  with  porosity.  In  Figure  13(c) 
it  appears  that  U  is  a  better  estimate  for  urand  and  u£  at  9  =  0.98  than  it  was 
for  9  =  0.95  depicted  in  Figure  12(c)  and  9  =  0.9  depicted  in  Figure  11,  so 
again  we  see  that  U  is  a  better  approximation  at  lower  levels  of  porosity  or 
higher  values  of  9. 
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Random  Sample,  8  =  0.9 
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Figure  11:  An  example  simulation  at  6  =  0.9,  or  equivalently  10%  poros¬ 
ity  (a)  The  random  geometry  £4;  (b)  The  geometric  approximation  of  £4 
assumed  in  homogenization  theory,  £4e  where  e  =  0.14;  (c)  urand(x,  0,  t), 
ue(x,0,t),  and  U(x,Q,t)  for  t  =  0.1  sec  on  024.  For  a  better  repre¬ 
sentation  of  the  behavior  of  the  three  systems  over  time  on  the  source 
boundary  (oq),  see  the  corresponding  movie  by  clicking  here.. 
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Random  Sample,  9  =  0.95 
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Figure  12:  An  example  simulation  at  9  =  0.95,  or  equivalently  5%  poros¬ 
ity  (a)  The  random  geometry  0;  (b)  The  geometric  approximation  of  Q 
assumed  in  homogenization  theory,  where  £  =  0.17;  (c)  wrand(x,  0,  t), 
u£(x,0jt),  and  U(x,0,t)  for  t  =  0.1  sec  on  0J4.  For  a  better  repre¬ 
sentation  of  the  behavior  of  the  three  systems  over  time  on  the  source 
boundary  (oj 4),  see  the  corresponding  movie  by  clicking  here. 
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Random  Sample,  0  =  0.98 
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Figure  13:  An  example  simulation  at  9  =  0.98,  or  equivalently  2%  poros¬ 
ity  (a)  The  random  geometry  fl;  (b)  The  geometric  approximation  of  0 
assumed  in  homogenization  theory,  where  e  =  0.17;  (c)  urand(a:,  0,  t), 
u£(x,0,t),  and  U(x,0,t)  for  t  =  0.1  sec  on  u;4.  For  a  better  repre¬ 
sentation  of  the  behavior  of  the  three  systems  over  time  on  the  source 
boundary  (oq),  see  the  corresponding  movie  by  clicking  here. 
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Model 

6  =  0.9 

9  =  0.95 

6  =  0.98 

Mean  time 

^rand 

45  sec 

53  sec 

62  sec 

u£ 

22  sec 

28  sec 

37  sec 

U 

3.9  sec 

4.0  sec 

3.7  sec 

Mean  mesh  size 

^rand 

7851 

8551 

9965 

u£ 

5016 

5744 

7317 

U 

355 

355 

355 

Table  4:  In  the  first  three  rows,  the  last  three  columns  contain  the 
average  time  to  calculate  the  respective  solution  for  the  200  simulations. 
The  last  three  rows  contain  the  average  mesh  size  over  200  simulations. 


We  again  carried  out  200  simulations  at  each  porosity  level  (8  =  0.9,  0.95 
and  0.98)  like  those  used  to  generate  Figures  11-13.  For  each  simulation,  we 
randomly  generated  the  domain  12  on  which  we  solved  (18)  for  urand.  We 
also  created  a  random  !2£,  to  solve  (20)  for  u£ .  We  used  the  reference  cell 
used  to  generate  12e  to  solve  (23)  for  U.  In  Table  4  we  report  the  average 
time  to  solve  each  of  the  systems  and  the  average  mesh  size  used  to  discretize 
the  respective  domain.  It  is  clear  that  of  the  three  systems  (18)  is  the  most 
computationally  intensive.  The  limit  system  (23)  is  the  least  computationally 
intensive;  its  mesh  is  an  order  of  magnitude  smaller  than  the  meshes  for  the 
other  two  systems.  Also  the  computing  time  required  to  solve  for  U  is  five  to 
ten  times  less  than  the  time  to  solve  for  u£  and  ten  to  fifteen  times  less  time 
to  solve  for  urand.  This  supports  the  use  of  the  limit  system  (23)  as  a  more 
computationally  suitable  model  than  (18)  and  (20). 

The  limit  system  (23)  seems  to  be  a  good  approximation  of  (18)  and  (20). 
We  consider  the  Frobenius  norm  of  the  difference  between  the  limit  systems 
on  the  source  boundary  oq.  We  will  take  values  at  time  points  tj  =  for 
j  G  {1, 2, . . .  101}  and  nodes  Xi  =  for  i  G  {1, 2, . . .  41}.  For  each  simulation, 
we  calculated  the  Frobenius  norm  of  \\U  —  urand\\p 

II u  -  urand||F  =  ^/e5  1  £5=1  {U{xi,  0,  tj)  -  u^d(Xi,  0,  tj))2,  (28) 

and  \\U  —  u£\\f-  We  report  the  average  and  standard  deviation  of  these  quan- 
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9  =  0.9 

9  =  0.95 

9  =  0.98 

Mean  |  wrand  —  U\\F 

1.21 

0.631 

0.267 

St.  Dev.  ||wrand  -  U\\F 

0.245 

0.109 

0.042 

Mean  | \u£  —  U\\F 

1.322 

0.656 

0.267 

St.  Dev.  | \u£  —  U\\f 

0.258 

0.109 

0.042 

Table  5:  The  first  and  third  row  contain  the  mean  Frobenius  norm  of 
the  difference  of  the  solutions  for  the  200  simiulations.  The  second  and 
fourth  row  contain  the  standard  deviation  of  the  Frobenius  norm  of  the 
difference  of  the  solutions  for  the  200  simulations. 


tities  for  each  porosity  level  in  Table  5.  As  8  increases,  U  better  approximates 
both  urand  and  u£ .  However,  it  is  not  clear  that  U  better  approximates  ei¬ 
ther  ue  or  urand  which  is  interesting  to  note  but  would  need  a  more  analytical 
treatment  to  investigate. 

5  Conclusion  and  Future  Work 

We  have  selected  and  analyzed  two  models,  which  we  can  use  to  represent 
a  flash- heat  experiment  mathematically.  The  first  model,  system  (9),  corre¬ 
sponds  to  a  small  loss  across  the  surface  of  the  cross  section,  and  the  second, 
system  (13),  corresponds  to  imperfect  insulation  leading  to  small  loss  on  all  of 
the  boundary  of  the  cross  section.  We  plan  to  continue  to  use  these  in  future 
investigations.  In  our  future  efforts  though,  we  will  focus  on  (13)  to  model  the 
flash-heat  experiment  as  it  is  a  more  physically  relevant  model  than  (9).  The 
results  of  homogenization  theory  in  Section  4  are  encouraging.  It  is  also  clear 
in  Section  4  that  the  choice  of  reference  domain,  specifically  the  location  of 
pores  in  fL  in  relation  to  the  source  boundary  W4  determines  the  efficacy  of 
U  to  approximate  us .  Similarly,  the  location  of  the  pores  in  Ll  in  the  random 
domain  determine  the  behavior  of  urand  in  relation  to  the  approximation  U. 
We  also  observed  that  U  approximates  urand  best  at  the  lowest  level  of  poros¬ 
ity,  98%  or  6  =  0.98.  In  the  future  work  we  would  like  to  use  limit  system 
(23)  developed  in  Section  3  as  a  model  to  carry  out  parameter  estimation  on 
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